//#ist=0x8088cac) at /usr/include/c++/4.6/bits/stl_construct.h:127
#include "CSpaceCharge.h"
#include <iostream>
#include <cmath>
#include <stdio.h>
#include <iomanip>
#include <sstream>
#include <iterator>
using namespace std;
const double dzero=1e-9;
const double dPermit=8.85e-12;
const double dElecCharge=1.602e-19;
int xtimes=1;
int ytimes=1;
//#include "common.h"

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])
//电流单位是mA
CSpaceCharge::CSpaceCharge()
{}
CSpaceCharge::~CSpaceCharge()
{
  delete dEnergyK;
  /*
     for(int iTmp=0;iTmp<nr;iTmp++)
     {
     delete dDensity[iTmp];
     delete dDensityK[iTmp];
     }
     delete dDensity;
     delete dDensityK;
   */
}
void CSpaceCharge::formGrid()
{
  /*
     plim(&dPhis, &dpLeft);
     double pright=-dPhis;
     double cay=2*pi/dLCell;
     double zmax=(pright-dpLeft)/cay;  
     double ddz=zmax/nz;
  //    if(dz<ddz)dz=ddz;    //dz is the long grid length,while ddz is the value after regulate the grid
   */
  vvvdaa.resize(nx);
  vvvdex.resize(nx);
  vvvdey.resize(nx);
  vvvdez.resize(nx);
  compData.resize(nx);
  for(int i=0;i<vvvdaa.size();i++)
  {
    vvvdaa[i].resize(ny);
    vvvdex[i].resize(ny);
    vvvdey[i].resize(ny);
    vvvdez[i].resize(ny);
    compData[i].resize(ny);
    for(int j=0;j<vvvdaa[i].size();j++)
    {
      vvvdaa[i][j].resize(nz);
      vvvdex[i][j].resize(nz);
      vvvdey[i][j].resize(nz);
      vvvdez[i][j].resize(nz);
      compData[i][j].resize(2*nz);
      for(int k=0;k<vvvdaa[i][j].size();k++)
      {
        vvvdaa[i][j][k]=0;
        vvvdex[i][j][k]=0;
        vvvdey[i][j][k]=0;
        vvvdez[i][j][k]=0;
        compData[i][j][2*k]=0;
        compData[i][j][2*k+1]=0;

      }
    }  
  }

  vdxm.resize(nx-1);
  vdym.resize(ny-1);
  vdzm.resize(nz-1);
  vdxs.resize(nx);
  vdys.resize(ny);
  vdzs.resize(nz);
  xHalf=(nx-1)*dx/2.0;
  yHalf=(ny-1)*dy/2.0;
  zHalf=(nz-1)*dz/2.0;
  for(int i=0;i<nx-1;i++)
  {
    vdxs[i]=i*dx;
    vdxm[i]=vdxs[i]+dx/2.0;
  }
  vdxs[nx-1]=dx*(nx-1);
  for(int i=0;i<ny-1;i++)
  {
    vdys[i]=i*dy;
    vdym[i]=vdys[i]+dy/2.0;
  }
  vdys[ny-1]=dy*(ny-1);
  for(int i=0;i<nz-1;i++)
  {
    vdzs[i]=i*dz;
    vdzm[i]=vdzs[i]+dz/2.0;
  }
  vdzs[nz-1]=dz*(nz-1);
}
/*
   rmax=nr*dr;
   pl =dBetaS*dWaveLength;
   vdrm.resize(nr+1); 
   vdrs.resize(nr);	
   vdrssq.resize(nr); 
   vdzm.resize(nz+1);
   vdzzs.resize(nz+1); 
   vders.resize(nr*(1+nz)*(nr+1));   //z处nr个源点在nz*(1+nr)个场点的电场
   vdezs.resize(nr*(1+nz)*(nr+1));   //
//	vder.resize((nr+1)*(nz+1));   //每个场点网格点的电场
//	vdez.resize((nr+1)*(nz+1));
vvdaa.resize(nr);
vvder.resize(nr+1);
vvdez.resize(nr+1);
for(int i=0;i<vvdaa.size();i++)
vvdaa[i].resize(nz+1);
for(int i=0;i<vvder.size();i++)
{
vvder[i].resize(nz+1);
vvdez[i].resize(nz+1);
}
vdrm[0]=0;
for (int i=1; i<nr+1; i++)
{
if(iflag==0)
{
//rm = 径向nr层网格每层的半径。从中心开始，rm[0]=0，直到最外层。
vdrm[i]=i*dr;
vdrssq[i-1]=0.5*(pow(vdrm[i-1],2)+pow(vdrm[i],2));
vdrs[i-1]=sqrt(vdrssq[i-1]);
}
else if(iflag==1)
{
vdrm[i]=i*dr-rHalf;
vdrs[i-1]=0.5*(vdrm[i-1]+vdrm[i]);
vdrssq[i-1]=vdrs[i-1]*vdrs[i-1];
}
//			rssq[nr]
}               
for(int i=0; i<nz+1; i++)
{
vdzm[i]=i*dz;
//			zm = 纵向每层网格的纵向坐标。zm[0]=0
vdzzs[i]=vdzm[i]+dz/2;
}
}
 */
void CSpaceCharge::scheff(double dist)
{
  //产生mesh
  dt=dist;
  dQ =dI/(dFreq*iParticleNumTotal);  // I is average current in one rf period
  double xbar=0;	//MainCode.h
  double ybar=0;	//MainCode.h
  double xsq=0;
  double ysq=0;
  double pbar=0;	//pbar 为正离子纵向坐标z的平均值	//MainCode.h
  double pmax=0;
  double pmin=0;

  for(int np=0; np<cord->size(); np++)
  {
    xbar=xbar+(*cord)[np][0];
    ybar=ybar+(*cord)[np][2];
    xsq=xsq+(*cord)[np][0]*(*cord)[np][0];
    ysq=ysq+(*cord)[np][2]*(*cord)[np][2];
    pbar=pbar+(*cord)[np][4];
    if((*cord)[np][4]<pmin) pmin=(*cord)[np][4];
    if((*cord)[np][4]>pmax) pmax=(*cord)[np][4];
  }  
  /*
     double plT=(pmax-pmin);
     double zlT=(pmax-pmin)/2/pi*dWaveLength*dBetaS;
     if(dz*nz>3*zlT) 
     {dpLeft=pmin-plT;dz=3*zlT/nz; 
     for(int i=0; i<nz+1; i++)
     {
     vdzm[i]=i*dz;
     vdzzs[i]=vdzm[i]+dz/2;
     }
     }
   */
  ng=cord->size();                                     //Rx---Right-----i
  xbar=xbar/ng;                                        //Lx---Left------i-1
  ybar=ybar/ng;                                        //Uy---Up--------j
  xsq=xsq/ng;                                          //Dy---Down------j-1
  ysq=ysq/ng;                                          //Fz---forward---k 
  pbar=pbar/ng;                                        //Bz---BackWard--k-1

  double dxP,dyP,dzP,dphiP;
  int    ix1,ix2,iy1,iy2,iz1,iz2;
  double dx1,dx2,dy1,dy2,dz1,dz2;
  for(int np=0;np<cord->size();np++)
  {
    dxP=(*cord)[np][0]+xHalf; dyP=(*cord)[np][2]+yHalf;dphiP=(*cord)[np][4];
    dzP=dphiP*dBetaS*dWaveLength/2.0/pi+zHalf;
    ix1=dxP/dx;iy1=dyP/dy;iz1=dzP/dz;
    ix2=ix1+1;iy2=iy1+1;iz2=iz1+1;

    if(ix1<0) {ix1=0;ix2=0;}
    if(ix1>=nx-1){ix1=nx-1;ix2=nx-1;}
    if(iy1<0) {iy1=0;iy2=0;}
    if(iy1>=ny-1){iy1=ny-1;iy2=ny-1;}
    if(iz1<0) {iz1=0;iz2=0;}
    if(iz1>=nz-1){iz1=nz-1;iz2=nz-1;}

    dx2=1.0;dy2=1.0;dz2=1.0;
    if(ix1!=ix2)
    {
      dx2=(dxP-vdxs[ix1])/(vdxs[ix2]-vdxs[ix1]);
    }
    if(iy1!=iy2)
    {
      dy2=(dyP-vdys[iy1])/(vdys[iy2]-vdys[iy1]);
    }
    if(iz1!=iz2)
    {
      dz2=(dzP-vdzs[iz1])/(vdzs[iz2]-vdzs[iz1]);
    }
    dx1=1.0-dx2;
    dy1=1.0-dy2;
    dz1=1.0-dz2;
    vvvdaa[ix1][iy1][iz1]=vvvdaa[ix1][iy1][iz1]+dx1*dy1*dz1;
    vvvdaa[ix1][iy1][iz2]=vvvdaa[ix1][iy1][iz2]+dx1*dy1*dz2;
    vvvdaa[ix1][iy2][iz1]=vvvdaa[ix1][iy2][iz1]+dx1*dy2*dz1;
    vvvdaa[ix1][iy2][iz2]=vvvdaa[ix1][iy2][iz2]+dx1*dy2*dz2;
    vvvdaa[ix2][iy1][iz1]=vvvdaa[ix2][iy1][iz1]+dx2*dy1*dz1;
    vvvdaa[ix2][iy1][iz2]=vvvdaa[ix2][iy1][iz2]+dx2*dy1*dz2;
    vvvdaa[ix2][iy2][iz1]=vvvdaa[ix2][iy2][iz1]+dx2*dy2*dz1;
    vvvdaa[ix2][iy2][iz2]=vvvdaa[ix2][iy2][iz2]+dx2*dy2*dz2;
  }
  for(int i=0;i<vvvdaa.size();i++)
    for(int j=0;j<vvvdaa[i].size();j++)
      for(int k=0;k<vvvdaa[i][j].size();k++)
      {
        vvvdaa[i][j][k]=vvvdaa[i][j][k]*dQ/dx/dy/dz*1.0e6;
        compData[i][j][2*k]=vvvdaa[i][j][k];
        compData[i][j][2*k+1]=0;
      }
  EFieldMap();
  ofstream of1("potentialXY");
  for(int i=0;i<nx;i++)
    for(int j=0;j<ny;j++)
      of1<<vdxs[i]-xHalf<<" "<<vdys[j]-yHalf<<" "<<(compData[i][j][nz-2]+compData[i][j][nz])/2.0<<endl;
  of1.close();
  of1.open("potentialYZ");
  for(int i=0;i<ny;i++)
    for(int j=0;j<nz;j++)
      of1<<vdys[i]-yHalf<<" "<<vdzs[j]-zHalf<<" "<<(compData[nx/2][i][2*j]+compData[nx/2-1][i][2*j])/2.0<<endl;
  of1.close();
  of1.open("potentialXZ");
  for(int i=0;i<nx;i++)
    for(int j=0;j<nz;j++)
      of1<<vdxs[i]-xHalf<<" "<<vdzs[j]-zHalf<<" "<<(compData[i][ny/2][2*j]+compData[i][ny/2-1][2*j])/2.0<<endl;
  of1.close();
  FieldError();
  /*
     ofstream ofX("aaX");
     double dTmp=0;
     for(int i=0;i<vvvdaa.size();i++)
     {
     for(int j=0;j<vvvdaa[i].size();j++)
     for(int k=0;k<vvvdaa[i][j].size();k++)
     dTmp+=vvvdaa[i][j][k];
     ofX<<vdxs[i]<<" "<<dTmp<<endl;;
     dTmp=0;
     }
     ofX.close();
     ofstream ofY("aaY");
     for(int i=0;i<vvvdaa[0].size();i++)
     {
     for(int j=0;j<vvvdaa.size();j++)
     for(int k=0;k<vvvdaa[j][i].size();k++)
     dTmp+=vvvdaa[j][i][k];
     ofY<<vdys[i]<<" "<<dTmp<<endl;;
     dTmp=0;
     }
     ofY.close();
     ofstream ofZ("aaZ");
     for(int i=0;i<vvvdaa[0][0].size();i++)
     {
     for(int j=0;j<vvvdaa.size();j++)
     for(int k=0;k<vvvdaa[j].size();k++)
     dTmp+=vvvdaa[j][k][i];
     ofZ<<vdzs[i]<<" "<<dTmp<<endl;;
     dTmp=0;
     }
     ofZ.close();
   */
  //  getVvv(&vvvdaa,"vvvdaa",0);
  //  getVvv(&compData,"potential",1,dx,dy,dz);
  //  getVvv(&vvvdaa,"vvvdaaYZ",0,0,0,0,1);
  //  getVvv(&vvvdaa,"vvvdaaXY",0,0,0,0,2);
  //  getVvv(&vvvdaa,"vvvdaaXZ",0,0,0,0,3);
  ////  getVvv(&compData,"potentialYZ",1,dx,dy,dz,1);
  ////  getVvv(&compData,"potentialXY",1,dx,dy,dz,2);
  ////  getVvv(&compData,"potentialXZ",1,dx,dy,dz,3);
  //  getVvv(&vvvdex,"vvvdexYZ",0,dx,dy,dz,1);
  //  getVvv(&vvvdex,"vvvdexXY",0,dx,dy,dz,2);
  //  getVvv(&vvvdex,"vvvdexXZ",0,dx,dy,dz,3);
  //  getVvv(&vvvdey,"vvvdeyYZ",0,dx,dy,dz,1);
  //  getVvv(&vvvdey,"vvvdeyXY",0,dx,dy,dz,2);
  //  getVvv(&vvvdey,"vvvdeyXZ",0,dx,dy,dz,3);
  //  getVvv(&vvvdez,"vvvdezYZ",0,dx,dy,dz,1);
  //  getVvv(&vvvdez,"vvvdezXY",0,dx,dy,dz,2);
  //  getVvv(&vvvdez,"vvvdezXZ",0,dx,dy,dz,3);
}

void CSpaceCharge::EFieldMap()      //求解每个场点网格点电场值
{
  FFT3D(&compData);
  //getVvv(&compData,"comp-before-FFT-check",1,0,0,0,0);
  /*
     ofstream of("comp-before-FFT-check");
     for(int i=0;i<compData.size();i++)
     for(int j=0;j<compData[i].size();j++)
     for(int k=0;k<compData[i][j].size();k=k+2)
     of<<vdxs[i]-xHalf<<" "<<vdys[j]-yHalf<<" "<<vdzs[k/2]-zHalf<<" "<<compData[i][j][k]<<endl;
     of.close();
   */
  //FFT1D(&compData);
  double dAngle;
  for(int i=0;i<compData.size();i++)
  {
    for(int k=0;k<compData[i].size();k++)
      for(int j=0;j<compData[i][k].size();j=j+2)
      {
        dAngle=atan2(compData[i][k][j+1],compData[i][k][j]);
        compData[i][k][j]=sqrt(compData[i][k][j]*compData[i][k][j]+compData[i][k][j+1]*compData[i][k][j+1]);
        compData[i][k][j+1]=dAngle;
      }
  }
  Potential();              
  double dAmp;
  for(int i=0;i<compData.size();i++)
  {
    for(int k=0;k<compData[i].size();k++)
      for(int j=0;j<compData[i][k].size();j=j+2)
      {
        dAngle=compData[i][k][j+1];
        dAmp=compData[i][k][j];
        compData[i][k][j+1]=dAmp*sin(dAngle);
        compData[i][k][j]=dAmp*cos(dAngle);
      }
  }
  IFFT3D(&compData);
  //  IFFT1D(&compData);
  for(int i=0;i<compData.size();i++)
  {
    for(int k=0;k<compData[i].size();k++)
      for(int j=0;j<compData[i][k].size();j=j+2)
      {
        dAngle=atan2(compData[i][k][j+1],compData[i][k][j]);
        compData[i][k][j]=sqrt(compData[i][k][j]*compData[i][k][j]+compData[i][k][j+1]*compData[i][k][j+1]);
        compData[i][k][j+1]=dAngle;
      }
  }
  //getVvv(&compData,"comp-after-IFFT-check",1,0,0,0,0);
  ofstream ofPotential("potential.txt");
  for(int i=0;i<vdxs.size();i++)
    for(int j=0;j<vdys.size();j++)
      for(int k=0;k<vdzs.size();k++)
        ofPotential<<vdxs[i]-xHalf<<" "<<vdys[j]-yHalf<<" "<<vdzs[k]-zHalf<<"  "<<compData[i][j][2*k]<<endl;
  ofPotential.close();

  stringstream ss;
  string sTmp;
  string filename;
  ss<<nx;
  ss>>sTmp;
  filename+=sTmp;
  filename+="-";
  ss.clear();
  ss<<dx;  ss>>sTmp;
  filename+=sTmp;
  sTmp=filename;
  filename+="-pot.txt";
  ofPotential.open(filename.c_str());
  ofPotential<<"n*d="<<nx/2<<"*"<<dx<<"  pot "<<endl;
  for(int i=0;i<vdxs.size();i++)
 //   ofPotential<<setw(8)<<vdxs[i]-xHalf<<"  "<<(compData[i][ny/2][nz]+compData[i][ny/2+1][nz]+compData[i][ny/2][nz+1]+compData[i][ny/2+1][nz+1])/4.0<<endl;
 ofPotential<<setw(8)<<vdxs[i]-xHalf<<"  "<<compData[i][ny/2][nz]<<endl;
  ofPotential.close();

  double dBUL,dBUR,dBDL,dBDR,dFUL,dFUR,dFDL,dFDR;
  for(int i=0;i<vdxm.size();i++)
    for(int j=0;j<vdym.size();j++)
      for(int k=0;k<vdzm.size();k++)                           //  zyx
      {                                                        //  BUL--back-up-left---(k--j+1--i)
        dBUL=compData[i][j+1][2*k];                            //  BUR--back-up-right--(k--j+1--i)
        dBUR=compData[i+1][j+1][2*k];
        dBDL=compData[i][j][2*k];
        dBDR=compData[i+1][j][2*k];
        dFUL=compData[i][j+1][2*(1+k)];
        dFUR=compData[i+1][j+1][2*(1+k)];
        dFDL=compData[i][j][2*(k+1)];
        dFDR=compData[i+1][j][2*(1+k)];
        vvvdex[i][j][k]=-(dBUR+dBDR+dFUR+dFDR-dBUL-dBDL-dFUL-dFDL)/4.0/dx*100;
        vvvdey[i][j][k]=-(dBUL+dBUR+dFUL+dFUR-dBDL-dBDR-dFDL-dFDR)/4.0/dy*100;
        vvvdez[i][j][k]=-(dFUL+dFUR+dFDL+dFDR-dBUL-dBUR-dBDL-dBDR)/4.0/dz*100;
      }
  ofstream ofField("efield.txt");
  for(int i=0;i<vvvdex.size()-1;i++)
    for(int j=0;j<vvvdex[i].size()-1;j++)
      for(int k=0;k<vvvdex[i][j].size()-1;k++)
        ofField<<vdxm[i]-xHalf<<" "<<vdym[j]-yHalf<<" "<<vdzm[k]-zHalf<<" "<<vvvdex[i][j][k]
          <<" "<<vvvdey[i][j][k]<<" "<<vvvdez[i][j][k]<<" "<<endl;
  ofField.close();

  //  ofField.open("error.txt");
  //  ofField.close();
  filename=sTmp;
  filename=sTmp+"-field.txt";
  ofField.open(filename.c_str());
  ofField<<"n*d="<<nx/2<<"*"<<dx<<"    e "<<endl;
  for(int i=0;i<vdxm.size();i++)
  {
//    if(vvvdex[i][ny/2-1][nz/2-1]>=0)
//      ofField<<vdxm[i]-xHalf<<" "<<sqrt(pow(vvvdex[i][ny/2-1][nz/2-1],2)
//          +pow(vvvdey[i][ny/2-1][nz/2-1],2)+pow(vvvdez[i][ny/2-1][nz/2-1],2))<<endl;
//    else
//      ofField<<vdxm[i]-xHalf<<" "<<-sqrt(pow(vvvdex[i][ny/2-1][nz/2-1],2)
//          +pow(vvvdey[i][ny/2-1][nz/2-1],2)+pow(vvvdez[i][ny/2-1][nz/2-1],2))<<endl;
//
    if(vvvdex[i][ny/2][nz/2]>=0)
      ofField<<vdxm[i]-xHalf<<" "<<sqrt(pow(vvvdex[i][ny/2][nz/2],2)
          +pow(vvvdey[i][ny/2][nz/2],2)+pow(vvvdez[i][ny/2][nz/2],2))<<endl;
    else
      ofField<<vdxm[i]-xHalf<<" "<<-sqrt(pow(vvvdex[i][ny/2][nz/2],2)
          +pow(vvvdey[i][ny/2][nz/2],2)+pow(vvvdez[i][ny/2][nz/2],2))<<endl;
  }
  ofField.close();
  Interpret();
}
void CSpaceCharge::FieldError()
{
  //建立一个球坐标系，
  double dGuss=0;
  double dphi,dtheta,dR;//单位是deg
  int iphi,itheta;
  double dxRadius,dyRadius,dzRadius;
  int ix1,ix2,iy1,iy2,iz1,iz2;
  double dx1,dx2,dy1,dy2,dz1,dz2;
  double ex,ey,ez;
  iphi=180;
  itheta=90;
  dphi=2*pi/iphi;
  dtheta=pi/itheta;
  dR=0.2;
  ofstream error("error.txt",ios::app);
  for(int i=0;i<itheta;i++)
    for(int j=0;j<iphi;j++)
    {
      dxRadius=dR*sin(i*dtheta)*cos(j*dphi);
      dyRadius=dR*sin(i*dtheta)*sin(j*dphi);
      dzRadius=dR*cos(i*dtheta);
      ix1=(dxRadius+xHalf)/dx;
      iy1=(dyRadius+yHalf)/dy;
      iz1=(dzRadius+zHalf)/dz;
      ix2=ix1+1;
      iy2=iy1+1;
      iz2=iz1+1;
      if(vdxm[ix1]>dxRadius&&ix1!=0&&ix1!=nx-1)
        ix2=ix1-1;
      else if((vdxm[ix1]>dxRadius&&ix1==0)||(vdxm[ix1]>dxRadius&&ix1==nx-1))
        ix2=ix1;
      if(vdym[iy1]>dyRadius&&iy1!=0&&iy1!=ny-1)
        iy2=iy1-1;
      else if((vdym[iy1]>dyRadius&&iy1==0)||(vdym[iy1]>dyRadius&&iy1==ny-1))
        iy2=iy1;
      if(vdzm[iz1]>dzRadius&&iz1!=0&&iz1!=nz-1)
        iz2=iz1-1;
      else if((vdzm[iz1]>dzRadius&&iz1==0)||(vdzm[iz1]>dzRadius&&iz1==nz-1))
        iz2=iz1;
      dx1=1.0;dy1=1.0;dz1=1.0;
      if(ix1!=ix2)
        dx1=(dxRadius+xHalf-vdxm[ix1])/(vdxm[ix1]-vdxm[ix2]);
      if(iy1!=iy2)
        dy1=(dyRadius+yHalf-vdym[iy1])/(vdym[iy1]-vdym[iy2]);
      if(iz1!=iz2)
        dz1=(dzRadius+zHalf-vdzm[iz1])/(vdzm[iz1]-vdzm[iz2]);
      dx2=1.0-dx1;dy2=1.0-dy1;dz2=1.0-dz1;
      ex=dz2*(dy2*(dx2*vvvdex[ix1][iy1][iz1]+dx1*vvvdex[ix2][iy1][iz1])
          +dy1*(dx2*vvvdex[ix1][iy2][iz1]+dx1*vvvdex[ix2][iy2][iz1]))
        +dz1*(dy2*(dx2*vvvdex[ix1][iy1][iz2]+dx1*vvvdex[ix2][iy1][iz2])
            +dy1*(dx2*vvvdex[ix1][iy2][iz2]+dx1*vvvdex[ix2][iy2][iz2]));
      ey=dz2*(dy2*(dx2*vvvdey[ix1][iy1][iz1]+dx1*vvvdey[ix2][iy1][iz1])
          +dy1*(dx2*vvvdey[ix1][iy2][iz1]+dx1*vvvdey[ix2][iy2][iz1]))
        +dz1*(dy2*(dx2*vvvdey[ix1][iy1][iz2]+dx1*vvvdey[ix2][iy1][iz2])
            +dy1*(dx2*vvvdey[ix1][iy2][iz2]+dx1*vvvdey[ix2][iy2][iz2]));
      ez=dz2*(dy2*(dx2*vvvdez[ix1][iy1][iz1]+dx1*vvvdez[ix2][iy1][iz1])
          +dy1*(dx2*vvvdez[ix1][iy2][iz1]+dx1*vvvdez[ix2][iy2][iz1]))
        +dz1*(dy2*(dx2*vvvdez[ix1][iy1][iz2]+dx1*vvvdez[ix2][iy1][iz2])
            +dy1*(dx2*vvvdez[ix1][iy2][iz2]+dx1*vvvdez[ix2][iy2][iz2]));
      dGuss+=dR*dR*1.0e-4*sin(i*dtheta)*dphi*dtheta*sqrt(ex*ex+ey*ey+ez*ez);
      //    error<<dGuss<<endl;
    }
  error<<dx<<" "<<dy<<" "<<dz<<" "<<nx<<" "<<ny<<" "<<nz<<" "<<(cord->size())
    <<" "<<dR<<" "<<dGuss<<" "<<1-dGuss/10.46<<endl;
  error.close();
}
void CSpaceCharge::FFT3D(vector<vector<vector<double> > > *vvvdComp)
{
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double aTmp[2*nz];
  for(int i=0;i<nx;i++)
  {              
    for(int j=0;j<ny;j++)
    {
      for(int k=0;k<nz;k++)
      {
        REAL(aTmp,k)=(*vvvdComp)[i][j][2*k];
        IMAG(aTmp,k)=(*vvvdComp)[i][j][2*k+1];
      }
      gsl_fft_complex_forward(aTmp,1,nz,compwave,compwork);
      for(int k=0;k<nz;k++)
      {
        (*vvvdComp)[i][j][2*k]=REAL(aTmp,k);
        (*vvvdComp)[i][j][2*k+1]=IMAG(aTmp,k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
  //getVvv(&compData,"comp-FFT-nx",0,0,0,0,0);
  compwork=gsl_fft_complex_workspace_alloc(ytimes*ny);
  compwave=gsl_fft_complex_wavetable_alloc(ytimes*ny);
  double bTmp[2*ny*ytimes];
  for(int i=0;i<nz;i++)
  {
    for(int j=0;j<nx;j++)
    {           
      for(int k=0;k<ny;k++)
      {
        REAL(bTmp,k)=(*vvvdComp)[j][k][2*i];
        IMAG(bTmp,k)=(*vvvdComp)[j][k][2*i+1];
      }
      for(int k=2*ny;k<2*ny*ytimes;k++)
      bTmp[k]=0;
      gsl_fft_complex_forward(bTmp,1,ny*ytimes,compwave,compwork);
      for(int k=0;k<ny;k++)
      {
        (*vvvdComp)[j][k][2*i]=REAL(bTmp,ytimes*k);
        (*vvvdComp)[j][k][2*i+1]=IMAG(bTmp,ytimes*k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
  //getVvv(&compData,"comp-FFT-nx-ny",0,0,0,0,0);

  compwork=gsl_fft_complex_workspace_alloc(nx*xtimes);
  compwave=gsl_fft_complex_wavetable_alloc(nx*xtimes);
  double cTmp[2*nx*xtimes];
  for(int i=0;i<nz;i++)
  {
    for(int j=0;j<ny;j++)
    {           
      for(int k=0;k<nx;k++)
      {
        REAL(cTmp,k)=(*vvvdComp)[k][j][2*i];
        IMAG(cTmp,k)=(*vvvdComp)[k][j][2*i+1];
      }
      for(int k=2*nx;k<2*nx*xtimes;k++)
      cTmp[k]=0;
      gsl_fft_complex_forward(cTmp,1,xtimes*nx,compwave,compwork);
      for(int k=0;k<nx;k++)
      {
        (*vvvdComp)[k][j][2*i]=REAL(cTmp,xtimes*k);
        (*vvvdComp)[k][j][2*i+1]=IMAG(cTmp,xtimes*k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
  //getVvv(&compData,"comp-FFT-nx--ny-nz",0,0,0,0,0);
}
void CSpaceCharge::FFT1D(vector<vector<vector<double> > > *vvvdComp)
{
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double aTmp[2*nz];
  for(int i=0;i<nx;i++)
  {              
    for(int j=0;j<ny;j++)
    {
      for(int k=0;k<nz;k++)
      {
        REAL(aTmp,k)=(*vvvdComp)[i][j][2*k];
        IMAG(aTmp,k)=(*vvvdComp)[i][j][2*k+1];
      }
      gsl_fft_complex_forward(aTmp,1,nz,compwave,compwork);
      for(int k=0;k<nz;k++)
      {
        (*vvvdComp)[i][j][2*k]=REAL(aTmp,k);
        (*vvvdComp)[i][j][2*k+1]=IMAG(aTmp,k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
}
void CSpaceCharge::IFFT3D(vector<vector<vector<double> > >* vvvdComp)
{
  //getVvv(&compData,"comp-before-IFFT",0,0,0,0,0);
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork=gsl_fft_complex_workspace_alloc(nx);
  compwave=gsl_fft_complex_wavetable_alloc(nx);
//  double aTmp[2*nx*xtimes];
  double aTmp[2*nx];
//  for(int i=0;i<2*nx*xtimes;i++)
//  aTmp[i]=0;
  for(int j=0;j<ny;j++)
  {
    for(int i=0;i<nz;i++)
    {
      for(int k=0;k<nx;k++)
      {
//        REAL(aTmp,xtimes*k)=(*vvvdComp)[k][j][2*i];
//        IMAG(aTmp,xtimes*k)=(*vvvdComp)[k][j][2*i+1];
        REAL(aTmp,k)=(*vvvdComp)[k][j][2*i];
        IMAG(aTmp,k)=(*vvvdComp)[k][j][2*i+1];
      }
      gsl_fft_complex_inverse(aTmp,1,nx,compwave,compwork);
      for(int k=0;k<nx;k++)
      {
        (*vvvdComp)[k][j][2*i]=xtimes*REAL(aTmp,k);
        (*vvvdComp)[k][j][2*i+1]=xtimes*IMAG(aTmp,k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
  //getVvv(&compData,"comp-IFFT-nx",0,0,0,0,0);
  compwork = gsl_fft_complex_workspace_alloc (ny);
  compwave = gsl_fft_complex_wavetable_alloc (ny);
  double bTmp[2*ny];
  //memset(bTmp,0,2*ny*ytimes*sizeof(double));
  for(int i=0;i<2*ny;i++)
  bTmp[i]=0;
  for(int i=0;i<nx;i++)
  {
    for(int j=0;j<nz;j++)
    {
      for(int k=0;k<ny;k=k++)
      {
        REAL(bTmp,k)=(*vvvdComp)[i][k][2*j];
        IMAG(bTmp,k)=(*vvvdComp)[i][k][2*j+1];
      }
      gsl_fft_complex_inverse(bTmp,1,ny,compwave,compwork);
      for(int k=0;k<ny;k++)
      {
        (*vvvdComp)[i][k][2*j]=ytimes*REAL(bTmp,k);
        (*vvvdComp)[i][k][2*j+1]=ytimes*IMAG(bTmp,k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
  //getVvv(&compData,"comp-IFFT-nx-ny",0,0,0,0,0);

  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double cTmp[2*nz];
  for(int i=0;i<nx;i++)
  {
    for(int j=0;j<ny;j++)
    {
      for(int k=0;k<nz;k=k++)
      {
        REAL(cTmp,k)=(*vvvdComp)[i][j][2*k];
        IMAG(cTmp,k)=(*vvvdComp)[i][j][2*k+1];
      }
      gsl_fft_complex_inverse(cTmp,1,nz,compwave,compwork);
      for(int k=0;k<nz;k++)
      {
        (*vvvdComp)[i][j][2*k]=REAL(cTmp,k);
        (*vvvdComp)[i][j][2*k+1]=IMAG(cTmp,k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
  gsl_fft_complex_wavetable_free(compwave);
  //getVvv(&compData,"comp-IFFT-nx-ny-nz",0,0,0,0,0);
}
void CSpaceCharge::IFFT1D(vector<vector<vector<double> > >* vvvdComp)
{
  gsl_fft_complex_wavetable *compwave;
  gsl_fft_complex_workspace *compwork;
  compwork = gsl_fft_complex_workspace_alloc (nz);
  compwave = gsl_fft_complex_wavetable_alloc (nz);
  double cTmp[2*nz];
  for(int j=0;j<ny;j++)
  {
    for(int i=0;i<nx;i++)
    {
      for(int k=0;k<nz;k=k++)
      {
        REAL(cTmp,k)=(*vvvdComp)[i][j][2*k];
        IMAG(cTmp,k)=(*vvvdComp)[i][j][2*k+1];
      }
      gsl_fft_complex_inverse(cTmp,1,nz,compwave,compwork);
      for(int k=0;k<nz;k++)
      {
        (*vvvdComp)[i][j][2*k]=REAL(cTmp,k);
        (*vvvdComp)[i][j][2*k+1]=IMAG(cTmp,k);
      }
    }
  }
  gsl_fft_complex_workspace_free(compwork);
}
void CSpaceCharge::Potential()       
{     
  double dkx[nx],dky[ny],dkz[nz],dif3K[nx][ny][nz];
  double dTmpx,dTmpy,dTmpz;
  double dif;
  dkx[0]=0;dky[0]=0;dkz[0]=0;

  //FFT3D======================================
  double dLz=dz*(nz-1);
  double dLx=dx*(nx-1);
  double dLy=dy*(ny-1);

  for(int k=0;k<ny;k++)
  {
    if(k<ny/2)
      dky[k]=k*(2*pi/dLy);
    else dky[k]=(k-ny)*(2*pi/dLy);
  }
  for(int i=0;i<nx;i++)           //when i=0.............!!!!!!!!
  {
    if(i<nx/2)
      dkx[i]=i*(2*pi/dLx);
    else dkx[i]=(i-nx)*(2*pi/dLx);
  }
  for(int j=0;j<nz;j++)
  {
    if(j<nz/2)
      dkz[j]=j*(2*pi/dLz);
    else       dkz[j]=(j-nz)*(2*pi/dLz);
  }
  /*
     vector<double> pulse;
     double dpulse
     itervvvdBegin=compData.begin();
     itervvvdEnd=compData.end();
     for(itervvvdBegin;itervvvdBegin!=itervvvdEnd;itervvvdBegin++)
     {
     itervvdBegin=itervvvdBegin->begin();
     itervvdEnd=itervvvdBegin->end();
     for(itervvdBegin;itervvdBegin!=itervvdEnd;itervvdBegin++)
     itervdBegin=itervvdBegin->begin();
     itervdEnd=itervvdBegin->end();
     for(itervdBegin;itervdBegin!=itervdEnd;itervdBegin++)
     {
     dTmp
     }
     }
   */
   /*
  vector<vector<vector<double> > >::iterator itervvvd;
  vector<vector<double> >::iterator itervvd;
  vector<double>::iterator itervd;

  itervvvdBegin=compData.begin();
  itervvvdEnd=compData.end();
  itervvvd=itervvvdBegin;
  int i,j,k;
  for(itervvvdBegin;itervvvdBegin!=itervvvdEnd;itervvvdBegin++)
  {
    i=0;
    itervvdBegin=itervvvdBegin->begin();
    itervvdEnd=itervvvdBegin->end();
    itervvd=itervvdBegin;
    for(itervvdBegin;itervvdBegin!=itervvdEnd;itervvdBegin++)
    {
      k=0;
      itervdBegin=itervvdBegin->begin();
      itervdEnd=itervvdBegin->end();
      itervd=itervdBegin;
      for(itervdBegin;itervdBegin!=itervdEnd;itervdBegin++)
      {
        j=0;
        if((itervvvdBegin==itervvvd)&&(itervvdBegin==itervvd)&&(itervdBegin==itervd))
          *itervdBegin=0;
        else
        {
          if(itervd==itervdBegin) dTmpz=1;
          else dTmpz=sin(dkz[j]*dz/2.0)/(dkz[j]*dz/2.0);
          if(itervvd==itervvdBegin) dTmpy=1;
          else dTmpy=sin(dky[k]*dy/2.0)/(dky[k]*dy/2.0);
          if(itervvvd==itervvvdBegin) dTmpx=1;
          else dTmpx=sin(dkx[i]*dx/2.0)/(dkx[i]*dx/2.0);
        }
        *itervdBegin=*itervdBegin/dif/dPermit;
      }
    }
  }
  */
  /*
  double dSumx,dSumy,dSumz;
  dSumy=0;dSumx=0;dSumz=0;
  for(int k=1;k<ny;k++)
  {
    dTmpy=sin(dky[k]*dy/2.0)/(dky[k]*dy/2.0);
    for(int i=1;i<nx;i++)           //when i=0.............!!!!!!!!
    {
      dTmpx=sin(dkx[i]*dx/2.0)/(dkx[i]*dx/2.0);
      for(int j=1;j<nz;j++)
      {
        dTmpz=sin(dkz[j]*dz/2.0)/(dkz[j]*dz/2.0);
        dif=dTmpx*dTmpx*dkx[i]*dkx[i]*1.0e4+dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4+dTmpy*dTmpy*dky[k]*dky[k]*1.0e4;
        compData[i][k][2*j]=compData[i][k][2*j]/dif/dPermit;
      }
    }
  }

//====test1,2D=====================
    for(int j=1;j<ny;j++)
    {
      for(int k=1;k<nz;k++)
        dSumz+=compData[0][j][2*k];
      compData[0][j][0]=-dSumz;
      dSumz=0;
    }
    for(int j=1;j<nz;j++)
    {
      for(int k=1;k<ny;k++)
        dSumy+=compData[0][k][2*j];
      compData[0][0][2*j]=-dSumy;
      dSumy=0;
    }

    for(int j=1;j<nx;j++)
    {
      for(int k=1;k<ny;k++)
        dSumy+=compData[j][k][0];
        compData[j][0][0]=-dSumy;
        dSumy=0;
    }
    for(int j=1;j<ny;j++)
    {
      for(int k=1;k<nx;k++)
        dSumy+=compData[k][j][0];
        compData[0][j][0]=-dSumy;
        dSumy=0;
    }

    for(int j=1;j<nx;j++)
    {
      for(int k=1;k<nz;k++)
        dSumy+=compData[j][0][2*k];
        compData[j][0][0]=-dSumy;
        dSumy=0;
    }
    for(int j=1;j<nz;j++)
    {
      for(int k=1;k<nx;k++)
        dSumy+=compData[k][0][2*j];
        compData[0][0][2*j]=-dSumy;
        dSumy=0;
    }

    for(int i=0;i<nx;i++)
    {
      for(int j=0;j<ny;j++)
      {
        for(int k=0;k<nz;k++)
          dSumy+=compData[i][0][2*k];
          dSumy+=compData[i][j][0];
          }
      compData[i][0][0]=-dSumy;
      dSumy=0;
    }

    for(int i=0;i<ny;i++)
    {
      for(int j=0;j<nx;j++)
      {
        for(int k=0;k<nz;k++)
          dSumy+=compData[0][i][2*k];
          dSumy+=compData[j][i][0];
      }
      compData[0][i][0]=-dSumy;
      dSumy=0;
    }


    for(int i=0;i<nz;i++)
    {
      for(int j=0;j<nx;j++)
      {
        for(int k=0;k<ny;k++)
          dSumy+=compData[0][k][2*i];
          dSumy+=compData[j][0][2*i];
      }
      compData[0][0][2*i]=-dSumy;
      dSumy=0;
    }

    for(int i=0;i<nz;i++)
      for(int j=0;j<nx;j++)
        for(int k=0;k<ny;k++)
          dSumy+=compData[j][k][2*i];
          compData[0][0][0]=-dSumy;


    */

  //====test2 -3D==================
  /*
  for(int i=1;i<nx;i++)
  {
    for(int j=1;j<ny;j++)
    {
      for(int k=1;k<nz;k++)
        dSumz+=compData[i][j][k];
      compData[i][j][0]=-dSumz;
      dSumz=0;
    }
    for(int j=1；j<nz;j++)
    {
      for(int k=0;k<ny;k++)
        dSumy+=compData[i][k][j];
    compData[i][0][2*j]=-dSumy;
    dSumy=0;
    }
  }
    */
  //=====orgion====================
      for(int k=0;k<ny;k++)
        for(int i=0;i<nx;i++)           //when i=0.............!!!!!!!!
          for(int j=0;j<nz;j++)
          {

            if((i==0)&&(j==0)&&(k==0))
              compData[i][k][2*j]=0;
            //compData[i][k][2*j]=0;
            else{
              if(k==0) dTmpy=1;
              else dTmpy=sin(dky[k]*dy/2.0)/(dky[k]*dy/2.0);
              if(i==0) dTmpx=1;
              else dTmpx=sin(dkx[i]*dx/2.0)/(dkx[i]*dx/2.0);
              if(j==0) dTmpz=1;
              else dTmpz=sin(dkz[j]*dz/2.0)/(dkz[j]*dz/2.0);
              dif=dTmpx*dTmpx*dkx[i]*dkx[i]*1.0e4+dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4+dTmpy*dTmpy*dky[k]*dky[k]*1.0e4;
              compData[i][k][2*j]=compData[i][k][2*j]/dif/dPermit;
            }
          }
    //====================================
    /*
    for(int k=1;k<ny;k++)
    {
      if(k<ny/2)
        dky[k]=k*(2*pi/dLy);
      else dky[k]=(k-ny)*(2*pi/dLy);
      dTmpy=sin(dky[k]*dy/2)/(dky[k]*dy/2);
      for(int i=1;i<nx;i++)           //when i=0.............!!!!!!!!
      {
        if(i<nx/2)
          dkx[i]=i*(2*pi/dLx);
        else dkx[i]=(i-nx)*(2*pi/dLx);
        dTmpx=sin(dkx[i]*dx/2)/(dkx[i]*dx/2);
        for(int j=1;j<nz;j++)
        {
          if(j<nz/2)
            dkz[j]=j*(2*pi/dLz);
          else       dkz[j]=(j-nz)*(2*pi/dLz);
          dTmpz=sin(dkz[j]*dz/2.0)/(dkz[j]*dz/2.0);
          dif=dTmpx*dTmpx*dkx[i]*dkx[i]*1.0e4+dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4+dTmpy*dTmpy*dky[k]*dky[k]*1.0e4;
          compData[i][k][2*j]=compData[i][k][2*j]/dif/dPermit;
        }                      
      }
    }
    */
    //getVvv(&compData,"comp-before-pot",0,0,0,0,0);
    //FFT2D============================================
    /*
       for(int i=0;i<nx;i++)
       for(int k=0;k<ny;k++)
       {
     if(k<ny/2)
     dky[k]=k*(2*pi/dy/ny);
     else dky[k]=(k-ny)*(2*pi/dy/ny);
     dTmpy=sin(dky[k]*dy)/(dky[k]*dy);

     for(int j=0;j<nz;j++)
     {
     if(j<nz/2)
     dkz[j]=j*(2*pi/dz/nz);
     else       dkz[j]=(j-nz)*(2*pi/dz/nz);
     dTmpz=sin(dkz[j]*dz)/(dkz[j]*dz);
     dif=dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4+dTmpy*dTmpy*dky[k]*dky[k]*1.0e4;
     compData[i][k][2*j]=compData[i][k][2*j]/dif/dPermit;
     }
     }
   */
  //FFT1D============================================
  /*
     double dL=dz*(nz-1);
     for(int j=0;j<nz;j++)
     {
     if(j<nz/2)
     dkz[j]=(j)*(2*pi/dL);
     else       dkz[j]=(j-nz)*(2*pi/dL);
     }
     for(int j=1;j<nz;j++)
     {
     dTmpz=sin(dkz[j]*dz/2.0)/(dkz[j]*dz/2.0);
     dif=dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4;
     for(int i=0;i<nx;i++)
     for(int k=0;k<ny;k++)
     compData[i][k][2*j]=compData[i][k][2*j]/dif/dPermit;
     }
     compData[0][0][0]=0;
   */
  //getVvv(&compData,"comp-after-pot",1,0,0,0,0);
}
void CSpaceCharge::setFreq(double _dFreq)
{
  dFreq=_dFreq*1e6;
  dWaveLength=clight/dFreq;
}
void CSpaceCharge::InputBeam(vector<vector<double> >* _vvdBeam,double _dI)
{
  cord=_vvdBeam;                                                    
  dI=_dI/1000.0;
  iParticleNumTotal=_vvdBeam->size();      //第一个粒子是同步粒子，不包括在粒子数范围内
  dEnergyK=new double[iParticleNumTotal];
  for(int n=0;n<iParticleNumTotal;n++)
    dEnergyK[n]=(*_vvdBeam)[n][5];
  double dEnergyS=dEnergyK[0];
  dGammaS=dEnergyS/dErest+1;
  dBetaS=sqrt(1-1/dGammaS/dGammaS);
  dLCell=dBetaS*dWaveLength;
}
void CSpaceCharge::InGrid(int _nx,int _ny,int _nz,double _dx,double _dy,double _dz)
{
  nx=_nx*2;           //in cm
  dx=_dx;
  ny=_ny*2;
  dy=_dy;
  nz=_nz*2;
  dz=_dz;
}
void CSpaceCharge::setPhis(double _dPhis)   //attention in degree
{
  dPhis=_dPhis/180*pi;
}
void CSpaceCharge::InputStruct(double _dFreq,double _dErest)
{
  dFreq=_dFreq*1.0e6;
  dWaveLength=clight/dFreq;
  dErest=_dErest;
}
void CSpaceCharge::Interpret()
{
  /*
     vector<vector<double> >::iterator iterBeginI=cord->begin();
     vector<vector<double> >::iterator iterEndI=cord->end();
     vector<double>::iterator iterBeginJ;
     vector<double>::iterator iterEndJ;
     for(iterBeginI;iterBeginI!=iterEndI;iterBeginI++)
     {
     iterBeginJ=iterBeginI->begin();
     iterEndJ=iterBeginI->end();
     }
   */
  double dxP,dyP,dzP;
  double dx1,dy1,dz1,dx2,dy2,dz2;
  int ix1,iy1,iz1,ix2,iy2,iz2;
  double dMass=dErest/clight/clight;
  double ex,ey,ez;
  for(int i=1;i<cord->size();i++)
  {
    dxP=(*cord)[i][0]+xHalf;
    dyP=(*cord)[i][2]+yHalf;
    dzP=(*cord)[i][4]*pi/180.0+zHalf;
    ix1=dxP/dx;
    iy1=dyP/dy;
    iz1=dzP/dz;
    ix2=ix1+1;
    iy2=iy1+1;
    iz2=iz1+1;
    if(vdxm[ix1]>dxP&&ix1!=0)
      ix2=ix1-1;
    else if(vdxm[ix1]>dxP&&ix1==0)
      ix2=ix1;
    if(vdym[iy1]>dyP&&iy1!=0)
      iy2=iy1-1;
    else if(vdym[iy1]>dyP&&iy1==0)
      iy2=iy1;
    if(vdzm[iz1]>dzP&&iz1!=0)
      iz2=iz1-1;
    else if(vdzm[iz1]>dzP&&iz1==0)
      iz2=iz1;
    dx1=1.0;dy1=1.0;dz1=1.0;
    if(ix1!=ix2)
      dx1=(dxP-vdxm[ix1])/(vdxm[ix2]-vdxm[ix1]);
    if(iy1!=iy2)
      dy1=(dyP-vdym[iy1])/(vdym[iy2]-vdym[iy1]);
    if(iz1!=iz2)
      dz1=(dzP-vdzm[iz1])/(vdzm[iz2]-vdzm[iz1]);
    dx2=1.0-dx1;dy2=1.0-dy1;dz2=1.0-dz1;
    ex=dz1*(dy1*(dx1*vvvdex[ix1][iy1][iz1]+dx2*vvvdex[ix2][iy1][iz1])
        +dy2*(dx1*vvvdex[ix1][iy2][iz1]+dx2*vvvdex[ix2][iy2][iz1]))
      +dz2*(dy1*(dx1*vvvdex[ix1][iy1][iz2]+dx2*vvvdex[ix2][iy1][iz2])
          +dy2*(dx1*vvvdex[ix1][iy2][iz2]+dx2*vvvdex[ix2][iy2][iz2]));
    ey=dz1*(dy1*(dx1*vvvdey[ix1][iy1][iz1]+dx2*vvvdey[ix2][iy1][iz1])
        +dy2*(dx1*vvvdey[ix1][iy2][iz1]+dx2*vvvdey[ix2][iy2][iz1]))
      +dz2*(dy1*(dx1*vvvdey[ix1][iy1][iz2]+dx2*vvvdey[ix2][iy1][iz2])
          +dy2*(dx1*vvvdey[ix1][iy2][iz2]+dx2*vvvdey[ix2][iy2][iz2]));
    ez=dz1*(dy1*(dx1*vvvdez[ix1][iy1][iz1]+dx2*vvvdez[ix2][iy1][iz1])
        +dy2*(dx1*vvvdez[ix1][iy2][iz1]+dx2*vvvdez[ix2][iy2][iz1]))
      +dz2*(dy1*(dx1*vvvdez[ix1][iy1][iz2]+dx2*vvvdez[ix2][iy1][iz2])
          +dy2*(dx1*vvvdez[ix1][iy2][iz2]+dx2*vvvdez[ix2][iy2][iz2]));
    //都是在靜止坐標系中考慮問題
    //相對粒子靜止的坐標系中，z=0，zp=0
    //假设beam.txt中数据存储的规则是，第一行是实验室坐标系同步粒子的坐标
    //其他行是束流坐标系的x、dx/dt,y,dy/dt,phi,dz/dt
    (*cord)[i][1]=ex*dQ*dt/dMass+(*cord)[i][1];
    (*cord)[i][3]=ey*dQ*dt/dMass+(*cord)[i][3];
    (*cord)[i][5]=dErest*dBetaS*pow(dGammaS,3)*ez*dQ*dt/dMass+(*cord)[i][5];
  }
  ofstream ofbeam("beam-step-3e-6");
  for(int i=0;i<cord->size();i++)
  {
  for(int j=0;j<(*cord)[i].size();j++)
  ofbeam<<(*cord)[i][j]<<" ";
  ofbeam<<endl;
  }
  // zpTozt(cord,dFreq);
}

void CSpaceCharge::getVvv(vector<vector<vector<double> > > *vvv,string str,int iflag,double dTmpx,double dTmpy,double dTmpz,int m)
  //compData iflag=1,else 0. dTmpx\dTmpy\dTmpz equals dx\dy\dz.  
  //if m=0, output all the datas in the vector,m=1 YZ plane;m=2 XY; m=3 XZ
{
  ofstream of(str.c_str());
  if(m==0)
  {
    for(int i=0;i<vvv->size();i++)
      for(int j=0;j<(*vvv)[i].size();j++)
      {
        if(iflag==0)  
          for(int k=0;k<(*vvv)[i][j].size();k++)
          {
            of<<vdxs[i]-xHalf+dTmpx/2.0<<" "<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[i][j][k]<<endl;
          }
        else if(iflag==1)
          for(int k=0;k<((*vvv)[i][j].size())/2;k++)
          {
            of<<vdxs[i]-xHalf+dTmpx/2.0<<" "<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[i][j][2*k]<<endl;
          }
      }
  }
  else if(m==1)   //输出y:z:vv[x==0][][]
  {
    for(int j=0;j<(*vvv)[nx/2-1].size();j++)
    {
      if(abs(vdxs[nx/2-1]-xHalf)<1e-5)
        //if(abs(vdxs[nx/2-1]-xHalf+dTmpx/2.0)<1e-5)
      {
        if(iflag==0)  
          for(int k=0;k<(*vvv)[nx/2-1][j].size();k++)
            of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[nx/2-1][j][k]<<endl;
        else if(iflag==1)
          for(int k=0;k<((*vvv)[nx/2-1][j].size())/2;k++)
            of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[nx/2-1][j][2*k]<<endl;
      }
      else 
      {
        if(iflag==0)  
          for(int k=0;k<(*vvv)[nx/2-1][j].size();k++)
            of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[nx/2-1][j][k]+(*vvv)[nx/2][j][k])/2.0<<endl;
        else if(iflag==1)
          for(int k=0;k<((*vvv)[nx/2-1][j].size())/2;k++)
            of<<vdys[j]-yHalf+dTmpy/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[nx/2-1][j][2*k]+(*vvv)[nx/2-1][j][k])/2.0<<endl;
      }
    }
  }
  else if(m==2)   //输出x:y:vv[x==0][][]
  {
    for(int j=0;j<(*vvv).size();j++)
    {
      if(abs(vdzs[nz/2-1]-zHalf)<1.0e-5||abs(vdzs[2*(nz/2-1)]-zHalf)<1.0e-5)     
        //if(abs(vdzs[nz/2-1]-zHalf+dTmpz/2.0)<1.0e-5||abs(vdzs[2*(nz/2-1)]-zHalf+dTmpz/2.0)<1.0e-5)     
      {
        if(iflag==0)  
          for(int k=0;k<(*vvv)[j].size();k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<(*vvv)[j][k][(nz/2-1)]<<endl;
        else if(iflag==1)
          for(int k=0;k<(*vvv)[j].size();k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<(*vvv)[j][k][2*(nz/2-1)]<<endl;
      }
      else 
      {
        if(iflag==0)  
          for(int k=0;k<(*vvv)[j].size();k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<((*vvv)[j][k][nz/2-1]+(*vvv)[j][k][nz/2])/2.0<<endl;
        else if(iflag==1)
          for(int k=0;k<(*vvv)[j].size();k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdys[k]-yHalf+dTmpy/2.0<<" "<<((*vvv)[j][k][2*(nz/2-1)]+(*vvv)[j][k][2*(nz/2)])/2.0<<endl;
      }
    }
  }
  else if(m==3)   //输出x:z:vv[x==0][][]
  {
    for(int j=0;j<(*vvv).size();j++)
    {
      if(abs(vdys[ny/2-1]-yHalf)<1.0e-5)
        //if(abs(vdys[ny/2-1]-yHalf+dTmpy/2.0)<1.0e-5)
      {
        if(iflag==0)  
          for(int k=0;k<(*vvv)[j][ny/2-1].size();k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[j][ny/2-1][k]<<endl;
        else if(iflag==1)
          for(int k=0;k<((*vvv)[j][ny/2-1].size())/2;k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<(*vvv)[j][ny/2-1][2*k]<<endl;
      }
      else 
      {
        if(iflag==0)  
          for(int k=0;k<(*vvv)[j][ny/2-1].size();k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[j][ny/2-1][k]+(*vvv)[j][ny/2][k])/2.0<<endl;
        else if(iflag==1)
          for(int k=0;k<((*vvv)[j][ny/2-1].size())/2;k++)
            of<<vdxs[j]-xHalf+dTmpx/2.0<<" "<<vdzs[k]-zHalf+dTmpz/2.0<<" "<<((*vvv)[j][ny/2-1][2*k]+(*vvv)[j][ny/2][2*k])/2.0<<endl;
      }
    }
  }
  of.close();
}
//void CSpaceCharge::error()
//{
//for(int i=0;i<vvvdex.size();i++)
//}
